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Abstract. We present particle-in-cell (PIC) simulations of minority energetic protons 
in deuterium plasmas, which demonstrate a collective instability responsible for 
emission near the lower hybrid frequency and its harmonics. The simulations capture 
the lower hybrid drift instability in a regime relevant to tokamak fusion plasmas, 
and show further that the excited electromagnetic fields collectively and collisionlessly 
couple free energy from the protons to directed electron motion. This results in an 
asymmetric tail antiparallcl to the magnetic field. We focus on obliquely propagating 
modes under conditions approximating the outer mid-plane edge in a large tokamak, 
through which there pass confined centrally born fusion products on banana orbits that 
have large radial excursions. A fully self-consistent electromagnetic relativistic PIC 
code representing all vector field quantities and particle velocities in three dimensions 
as functions of a single spatial dimension is used to model this situation, by evolving the 
initial antiparallcl travelling ring-beam distribution of 3MeV protons in a background 
lOkcV Maxwellian deuterium plasma with realistic ion-electron mass ratio. The 
simulations thus demonstrate a key building block of alpha channelling scenarios for 
burning fusion plasmas in tokamaks. 
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1. Introduction 

There is sustained interest in the collective instabilities of energetic ion populations in 
plasmas with confining magnetic fields that involve emission and absorption of radio 
frequency (RF) waves. For fusion plasmas, the primary aspects of practical interest 
include the suprathermal ion cyclotron emission arising from the products of deuterium 
and tritium fusion reactions, observed in both JET and TFTR, typically in the >10 
MHz range [T]|; and, at the highest frequencies where electron dynamics also enter, 
alpha channelling j2]. This refers to the possibility of collectively coupling energetic ion 
excited RF wave energy to create useful distortions of particle velocity distributions and 
hence, for example, to efficiently drive electron plasma current. In the present paper 
we focus on the second of these topics. As noted in the pioneering theoretical work of 
Fisch and co-workers [2] , understanding and exploiting the physics of lower hybrid waves 
that are excited collectively by the energetic ions and subsequently absorbed by Landau 
damping on the electrons, is central to several aspects of alpha channelling. This raises 
interesting theoretical challenges: in particular, the need to address simultaneously and 
self consistently the collective dynamics of ions, for which the highest natural frequency 
is the ion cyclotron frequency Q C i (and its harmonics), and the dynamics of waves 
involving collective motion of the electrons, for which the lowest natural frequency is 
the lower hybrid frequency pj], ujlh-, defined by 

i + n ce n a /io 2 p J 1 j 

Here Q ce , VL ci and u P i represent the electron and ion cyclotron frequencies and the ion 
plasma frequency respectively. Fusion alpha-particles are born at energies of 3.5MeV 
in deuterium-tritium fusion plasmas whose thermal energies are of order lOkeV. This 
creates the possibility of population inversion in velocity space, at least transiently 
during burn initiation before the alpha-particles have slowed through collisions with 
electrons. Population inversion opens the door to rapid energy exchange through 
collective instability. It is already known that, in tokamak plasmas, spatially localised 
inversions of the energy distribution of fusion-born ions can be sustained for long periods 
as a result of the interplay between particle energy and the pitch angle-dependent 
character of particle drift orbits. Specifically, ion cyclotron emission was excited in the 
outer mid-plane edge region of JET and TFTR as a consequence of the distinctive radial 
excursions (4j[5] of the drift orbits of fusion products born with pitch angles just inside the 
trapped-passing boundary. This was first seen for fusion protons born in pure deuterium 
plasmas in JET [4l|6], and subsequently for fusion alpha-particles in deuterium-tritium 



plasmas in JET [7j|8] and TFTR |1|[5||9], as well as injected beam ions 10 in TFTR and 
heated minority ions in JET [8]. For reviews and further theory, we refer to (Til, 12 . 
These observations motivate the choice of model distribution for the energetic ions in 
the present study. Given the variety of possible alpha channelling physics scenarios, it 
is necessary to find a focus for the representation of population inversion in the present 
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study, so we confine our attention to the primary type of ion distribution that has so far 
been observed to excite RF waves by collective instability in large tokamak edge plasmas. 
Specifically, we follow [9 11 13 in adopting a simple model in which the distribution of 
energies perpendicular to the local magnetic field is narrow, and there is a parallel drift. 
We note that similar distributions, giving rise to similar collective instability leading to 



excitation of RF waves, have been observed in the terrestrial magnetosphere 14 , 15 



and may be operative in the interstellar medium 16 18 . Our choice of parameter 



values for the model energetic ion distributions, and for the majority thermal ion and 
electron populations that comprise the background plasma, is similarly motivated by 
the observations of ion cyclotron emission from large tokamak edge plasmas. 

In this paper, we use a fully kinetic particle-in-cell (PIC) code to examine the 
collective excitation of waves near the lower hybrid frequency by such energetic ion 
populations in tokamak edge plasma conditions. Electrons damp the waves produced 
by the instability and are thereby accelerated. This results in an asymmetric tail in the 
distribution of electron velocities in the direction parallel to the magnetic field, implying 
an electron current. The fully kinetic PIC simulation self-consistently evolves in time 
all field vector quantities and particle velocities in three dimensions as a function of 
one space coordinate for three populations — energetic minority ions, thermal majority 
ions, and thermal electrons represented with realistic mass ratios. The results included 
in this paper expand and elaborate upon the mechanism identified in Ref. 19 for the 
generation of a suprathermal tail in the distribution of electron velocities parallel to 
the magnetic field and consequently provides a building block for the construction of 
tokamak-relevant alpha channelling scenarios. 

The linear phase of the instability underlying the phenomena seen in these fully 
nonlinear numerical simulations is consistent with a lower hybrid drift instability 
(LHDI). To simplify and generalise somewhat, the wave excited is predominantly 
electrostatic and propagates nearly perpendicular to the magnetic field, with a parallel 
phase velocity which (in the cases we consider) is sufficiently close to the electron thermal 
velocity to enable significant electron Landau damping. The wave frequency is high in 
comparison to the ion cyclotron frequency, so that the action of the energetic ions on 
the timescale of relevance is predominantly that of a quasi-unmagnetised beam. 

The scenario we consider differs fundamentally from that considered in much of the 
LHDI literature in one important aspect. Most previous work on LHDI has considered 
a free energy source residing in bulk ion drifts that are consistent with, and contribute 
to, the overall equilibrium of the plasma; notably, diamagnetic drifts. For the tokamak 
application, however, it is necessary to consider the LHDI of diffuse minority energetic 
ion populations that do not contribute to the plasma equilibrium. This is the case 
both for the edge localised population that appears to drive ion cyclotron emission, 
and for many alpha channelling scenarios. This degree of freedom arises from the 
interplay between particle energy and particle orbits in tokamak magnetic geometry, 
which ultimately derives from the combination of the particle toroidal velocity with 
the toroidal component of the magnetic vector potential in the canonical toroidal 
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momentum. It can enable the local spatial generation of anomalous velocity space 
structures. 

The LHDI is a phenomenon arising from drifting populations of ions, and its 
character is highly dependent on the values of key plasma parameters 20 23 . LHDI 



is present in natural 24 -26 and laboratory plasmas 27, 28 



LHDI with reference to magnetic reconnection in both laboratory plasmas 28 



Work has been done on 
and 



magnetospheric 29 31 plasmas. Inhomogeneities drive particle drifts in plasmas which 



generate the LHDI in the cases considered in most of the literature 23, 27 28 30 33 



whereas the much larger drift relative to the ion thermal velocity considered here is a 
consequence of the energy released per particle in fusion reactions and of the geometry 
of the confining magnetic field in the tokamak. A significant number of independent 
plasma parameters must be specified in calculations of the properties of the LHDI for 
a particular application. The quantitative characterisation of the LHDI in the research 
literature is thus inevitably incomplete, and will likely remain so, because it is possible 
to explore only patchily the multi-dimensional parameter space. This situation also 
applies as regards the nonlinear development of the LHDI. Davidson et al 23 consider 
the LHDI using imposed drifts from field inhomogeneities for the strictly perpendicular 
case. Hsia et al 22 includes finite k» in their analytical approach. The more recent 
study by Silveira et al 32 builds on these studies, and other important works too 
numerous to mention by name here. They present a unified formulation for instabilities 
in the lower hybrid range of frequencies due to field inhomogeneities. In these cases 
the LHDI nevertheless is excited by inhomogeneity induced drifts that are slower, in 
dimensionless terms, than the fusion product drifts considered in this paper. The 
pioneering computational study by Chen and Birdsall [34] has cold fluid electrons. 



Further numerical work by Yoon et al 21 investigates the LHDI with gyro-averaged 
electrons and kinetic ions, while work by Daughton et al 31 uses physical mass ratio 
kinetic particle simulations for Vn and V-B inhomogeneities which generate low drift 
velocities. 

Finally we note that PIC codes have been used to address astrophysical problems 
that involve the same combination of effects that is considered here: beam-type minority 
ion distributions, resulting in RF wave excitation that leads on, through damping, to 
electron acceleration 35 38 . PIC codes are also widely used to study other fundamental 



instabilities. Population inversion 39 
Kelvin-Helmholtz instability 



kink instability 40 , Weibel instability 41 



42 



shocks 43 44 and wave-particle interactions 45 -49 



are among the many avenues investigated using PIC. 



2. Simulations 



The 1D3V PIC code epochld, which is based on PSC [50], represents the full 
distributions of all the species in the simulation, using computational macroparticles 
in three dimensional velocity space and in one configuration space dimension. These are 
evolved in time via the relativistic Lorentz force law in the absence of collisions. From 
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this particle representation of the distribution function, the plasma moments (charge 
and current densities) are interpolated onto a fixed spatial grid in one dimension. From 
these, using the full Maxwell equations, the electromagnetic fields are advanced on the 
grid. The PIC algorithm thus self-consistently advances three dimensional vector fields 
and particle velocities in one spatial coordinate (x) and time. 

The spatial domain is split into No cells of width Ax, each of which has a grid point 
at the centre. Cell size, which is uniform across the domain, is set at 1/10 the electron 
Debye length and resolves the gyroradii of all species. The phase space density function 
/ of each species is represented by a large (682670) set of computational macroparticles 
which are initially distributed evenly in space amongst 2048 cells. Although all 
components of velocity are evolved, only the velocity component in the direction of the 
simulation domain causes the particles' position to change. Each velocity component 
of all particles contributes to the respective current component on the grid points in 
nearby cells. The field and particle advance timestep, At, divided by the cell size, Ax, 
resolves the speed of light, and thus resolves cyclotron frequencies. Periodic boundary 
conditions are used for particles and fields and any wavemodes present propagate parallel 
or antiparallel to the direction of the simulation domain. 

The trajectories of background thermal electrons and deuterons, both of which 
have non-drifting Maxwellian distributions with temperature of lOkeV, are evolved in 
time alongside the minority energetic protons (initially at 3MeV) which initially have a 
ring-beam velocity distribution modelled by 

f P = 2^ S ( V \\ -u)8(v± -v r ) (2) 

see for example Dendy et al. [13] and McClements et al. [9]. Here u is the velocity 
along the magnetic field and v r is the perpendicular velocity of the ring. As mentioned 
previously, this choice of model is motivated by multiple studies of ion cyclotron emission 
in JET and TFTR. For this reason the initial pitch angle of the ring distribution (angle 
between u and v r ) is initialised at 135° with respect to the background magnetic field 
at the start of the simulation. The initial velocity distributions of all the species in the 
simulation are uniform in space. The angle of inclination of the direction of spatial 
variation x, of our simulation to the initial background magnetic field is 9 = 84°. 
We simulate a quasi-neutral plasma with an electron number density n e = 10 18 m~ 3 , 
and an applied magnetic field of B — 3T. The electron beta is /3 e = 0.03%. The 
ratio of the initial energetic ion speed v p (t = 0) = (u 2 + f^) 1//2 , to the Alfven speed 
Va, is v p (t = 0)/Va — 0.51. This set of model parameters is selected to enable 
us to address quasi-perpendicular propagating modes under conditions approximating 
the edge plasma of a large tokamak, subject to the necessary simplifications. It 
optimizes the grid size with respect to the gyroradii of the particles, so as to resolve 
the essential physics with reasonable computational resources. The gyroradius of the 
energetic protons, and the gyroradii at thermal speeds of electrons and deuterons, are 
resolved by 965, 1.8 and 111 grid cells respectively. The ratio of protons to deuterons 
np/nd = 10~ 2 . The 3MeV ring of protons is not replenished as the distribution function 
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evolves away from the initial configuration; there is no source term. These parameter 
values imply that the total energy of the energetic proton population is 1.7 times that 
of the background thermal deuterons and electrons combined. This ratio exceeds by 
a factor ~ 10 the value anticipated for various energetic ion populations in next step 
fusion plasmas, see for example Table 1 of S. Putvinski (Nucl. Fus. 38, 1275 (1998)). 
This value is necessary to drive the instability on an acceptable timescale given finite 
computational resources, and brings this study closer to fusion relevant regimes than 
much of the LHDI literature. 



3. Flow of energy between particles and fields 

In order to examine the instability produced by the free energy in the energetic proton 
population we first consider the time evolution of the electric Ee and magnetic Eb total 
field energies in the simulation: 

Ng Ng 

^ = f£^ 2 A* e B = ^-J2^- B ^ Ax ( 3 ) 

i=l ^° i=l 



where B is the applied magnetic field. Figure 1(a) shows that the electric field energy, 
Ee, (light blue trace) starts rising exponentially after a period of ~10tlh- The magnetic 
field energy, Eb (magenta trace), rises at approximately the same rate as the electric 
field, but contains nearly an order of magnitude less energy, implying that the waves 
produced by the instability are largely electrostatic. The total kinetic energy of species 
i is defined by 

£j = Wj y^U/fffe 2 - mfc 4 - uijC 2 ) (4) 
i=i 

where Wi is the weighting factor to map the energy of each computational macroparticle, 
j, to the energy of the number of physical particles it represents. Particle weight is 
defined as Wi = n 0ti L/Np^ where n ,i is the species number density, L is the length 
of the spatial domain, and Npi is the number of computational macroparticles that 



represent the species in the simulation. The essential features of figure 1(a) indicate a 
linear phase of field growth during 10 < t/rLH < 15, accompanied by a corresponding 
small decline in total proton energy. Field amplitudes remain approximately constant 
after t > 20tlh, while proton energy continues to decline and is matched by an increase 
in electron energy. 

We now define a quantity which better shows the variations in kinetic energy of the 
proton population, which we denote as fluctuation energy e summed over the j = 1—Np^ 
computational macroparticles of each i th species: 

e,(f) <e,(f = 0)> |) (5) 
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The quantity e is the sum of the modulus of the change in kinetic energy of each particle 
from the initial rest frame mean energy value of the species, < £j(i = 0) >. Thus £j 
captures kinetic energy dispersion regardless of sign, whereas a time evolving change in 
total kinetic energy, AEi = Ei(t) — Ei{t = 0) where Ei(t = 0) = Yl!j=i £ j(t = 0)> om y 
represents changes in the combined energy of the ensemble of particles. In particular we 
use e to quantify the initial energy dispersion of the protons during the early phases of 
the instability, see figure l(b)| and the discussion in section |5j because in some respects 
it is less susceptible to noise. 

Figure 1(b) shows in greater detail the species energy dynamics. The values of 
proton fluctuation energy e p , and of the changes Ae e and Ae^ in total electron and 
deuterium kinetic energy are plotted alongside changes in ee and Eb- The value of e p 
grows, ultimately increasing by three orders of magnitude, whereas figure 1(a) shows 
that the total proton kinetic energy declines by much less than one order of magnitude. 
This reflects the role of the proton population as the source of free energy. The change 
in electron kinetic energy Ae c rises in close alignment with electric field energy ee 
during the linear stage of the simulation, 10 < t/rLH < 15. This is evidence of the 
electron oscillation within the largely electrostatic waves excited by the instability. The 
corresponding effect for the deuterons is also visible in figure l(b)| but is much less due 
to their higher mass. 




Figure 1. Panel (a): Time evolution of the kinetic energy, e,; denned by equation 
|1J in each plasma species, and of the energy of electric field Ee, and magnetic field 
Eb defined by equation [3j The energy of the magnetic field excludes the applied 
component. Panel (b): Time evolution of the proton fluctuation energy, ii defined by 
equation[5j change in kinetic energy of electrons and deuterons, Ae^ = £i(t) —£i(t = 0), 
and of the electric and magnetic field energies Ee and Eb- Both panels show vertical 
lines labelled (i) — (iv) which denote the four snapshots in time at which particle data 
is shown in other figures. Time is measured in units of the lower hybrid period. 
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Figure 2. Panel (a): Initial proton momentum PDF: along the simulation domain, p x \ 

in the plane created by the magnetic field and the simulation domain, p y ; perpendicular 

to both magnetic field and simulation domain p z = pj_ 2 ; parallel to the magnetic field 

<P\\\ and total perpendicular p± = (p j_ 1 +p1,2) 1 ^ 2 - Panel (b): Electron PDFs in 

relativistic parallel momentum , for the four snapshots in time (i) — (iv) defined in 

Fig. [I] also marked by vertical lines in panels (c) and (d). Vertical lines indicate 

the phase velocities of the forward and backward travelling dominant waves. Panel 

(c): Skewness 7 e of the electron vu PDF (black circles) and change in electron kinetic 

energy Ae e = e e (t) — e e (t = 0), (red crosses). Panel (d): Time evolving parallel 

current in the electron bulk v x > —2.5v e (t — 0) (red triangles), and in the electron tail 

v x < ~2.5v e (t — 0) (black stars). Momentum axes are in units of the initial relativistic 

T n - 1 1/2 

momentum pi(t = Q) = (l/Np t i)T,j_fpj.pj(t = 0) of the species i, and time axes 

are in units of the lower hybrid oscillation period. 



Figure 2(a) shows projections in momentum space of the PDF of the proton 
population at t = defined by (|2]). The ring-beam distribution, as aligned with respect 
to the magnetic field, is projected into the cartesian co-ordinates of the simulation where 
x is along the simulation domain, z is perpendicular to both the simulation domain and 
the magnetic field and y completes the orthogonal set. The initial PDF is a delta 
function in both pn and p± = [p\ 1 +£>j_ 2 ) 1 ^ 2 - The three momentum components are 
projected into the simulation domain, in which the ring distribution becomes a more 
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complicated function. Due to the oblique angle between simulation domain and applied 
magnetic field, the peaks in the PDF in p x (red circled trace) are not equal and opposite. 
This has consequences for the wave-particle interaction which are dealt with in the next 
section. 



Figure 2(b) shows electron parallel momentum at four different snapshots in time: 
(i) = 6.8tlh , (H) = 13.lT£ff, (Hi) = 16.8tx#, (iv) = 19. 9tlh- Henceforth, these 
snapshots in time will be referred to by their Roman numeral. 

Figures 2(b) and 2(c) show that for t > 15tlh the electron distribution develops 
an asymmetric tail in vn with corresponding skewness, 

\ 3 

vy(t)- < v u (t) > 



l(t) 



(6) 



where < v\\,i(t) > and a v .. . (t) are respectively the mean and standard deviation of electron 
v\\ at time t. 



Figure 2(d) shows the time evolving parallel electron current summed across 
configuration space. We plot the current for the bulk (v\\ > —2.5v e (t = 0)) and tail 
(v\\ < —2.5v e (t = 0)) of the electrons separately. We see that approximately equal 
and opposite currents are created by the bulk and tail electron populations in this 
collisionless simulation. This is a consequence of the periodic boundary conditions used 
in the code. The rate of change of current across the simulation domain under periodic 
boundary conditions follows from Ampere's and Faraday's laws: 



/'o 



d_ 
dt 



L Q2 

Vx Vx Edx - /j e — 



Jd.r = - I Vx Vx Ed.r - // (l r n — / Edx (7) 

where in one dimension V = (d/dx,0,0) and V x Vx = (0, — J^, — J^)- Periodic 
boundary conditions enforce F(0, t) = F(L, t) and V-F(0, t) = VF(L, t) where F is any 
component of the electromagnetic field. Equation is a wave equation in the plasma 
medium. For current density integrated over the box 



c 2 VxVx 



dt 2 



Jdx 



M 9 



Jdx 



(8) 



where J = crE and a is the plasma conductivity tensor. Consequently, any current 
density integrated across the domain may oscillate about zero but cannot grow in time 
for these simulations. We also note that collisions would differentially affect the bulk 
and tail electron populations, and hence their relative drift; however our PIC code does 
not incorporate collisions. 

The asymmetry in the electron parallel PDF, reflecting net directional electron 
acceleration, continues to grow in the period beyond t = 20tlh, during which wave 
energy quantified by Ee and Eb is approximately constant, see figure 



la 



We infer 



Landau damping of the excited waves on resonant electrons, which results in the 
flattening of the electron parallel momentum PDF in figure 2(b) The vertical lines in 
this panel indicate the electron nonrelativistic momentum that corresponds to the phase 
velocity of the dominant relevant features in the u, k transform of the electromagnetic 
field, to which we now turn. 
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4. The nature of the electromagnetic fields 

Figure [3] displays the spatio-temporal fast Fourier transform (FFT) of the electric 
field component along the simulation domain E x , for three epochs of the simulation, 
showing temporal evolution of the u, k structure. Referring to figure [TJ the three 
epochs correspond to the following phases of electromagnetic field evolution: initial, 
5 < t/r LH < 10; early linear growth, 10 < t/r LH < 15; and late linear growth, 
15 < t/rLH < 20. The peaks in intensity at k ~ ±2u pe /c and u ~ Qulh in all 
three panels of figure [3] reflect the resonant interaction of the energetic initial proton 
population with the normal mode of the background plasma to which they most easily 
couple, namely the electron-deuteron lower extraordinary- Whistler wave mode whose 
analytical characteristics are outlined below. 
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Figure 3. Spatio-temporal FFTs of the electric field component in the direction of 
the simulation domain E x , integrated over all x and the following time periods. Panel 
(a): 5 < t/r LH < 10. Panel (b): 10 < t/r LH < 15. Panel (c): 15 < t/r LH < 
20. Overplotted on all panels is the cold plasma dispersion relation (dash trace). 
Overplotted on panel (a) and panel (b) are velocities corresponding to the peaks in the 
initial proton v x PDF (dot-dash trace). Overplotted on panel (c) are the velocities of 
the peaks and extremes in the proton v x PDF at snapshot (iv) (dot-dash trace). The 
frequency axes are in units of lower hybrid frequency and wavenumber axes in inverse 
electron skin depths uj pe /c. 



Figure 3(a) shows the waves present in the simulation during its initial phase 
before the onset of linear growth of electric and magnetic fields, in the time interval 
5 < t/r LH < 10. Overplotted is the cold plasma electron-deuteron dispersion relation, 
implicitly defined by 

/ e 1 -N 2 -ie 2 \ 

is 2 e l -N 2 cos 2 {e) N 2 sin(6)cos(6) 
K N 2 sin(6)cos(6) E 3 -N 2 sm 2 (6) j 

This is the determinant of the cold plasma dielectric tensor for arbitrary propagation 
direction 9 with respect to the magnetic field, where the refractive index N = ck/u. 
The components of the dielectric tensor are 







(9) 
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where u P i and Q C i are the plasma and cyclotron frequencies of species i. 

The cold plasma dispersion relation shows predominantly perpendicular propagat- 
ing normal wavemodes; ordinary and extraordinary. However some properties of par- 
allel propagating waves (Whistler) are additionally present, because the angle of wave 
propagation is oblique. The result is an extraordinary- Whistler branch with a cut-off 
frequency > ujlh- The instability is still essentially a lower hybrid drift-type instability, 
although it is active at frequencies above ujlh- 



Figure 3(b) is obtained using data from the interval 10 < t/rtH < 15. Figure 



1(a) shows that this corresponds to the early linear stage of electric field wave growth, 
during which the changes in energy of electrons and deuterons also rise in tandem, 
reflecting their coherent participation in wave motion. The trace corresponding to 



snapshot (ii) in figure 2(b) shows that a monotonically decreasing electron tail is drawn 



out during this epoch, with skewness becoming apparent (figure 2(c)) at its end. Figure 



3(b) enables us to identify the u, k character of the waves excited in the linear growth 
phase of the simulation. Both forward and backward propagating waves are present, 



at Quilh % w ^ 8ulh and 2u pe /c < \k\ < 5u pe /c. Figure 3(b) also exhibits two 



peaks in electric field intensity which appear to be backward propagating harmonics 
of the fundamental activity at frequency u> ~ Qulh- During the late linear phase 
15 < t/r^H < 20, these multiple higher harmonics become established. 



Figure 3(c) shows data from 15 < t/rLH < 20 which encompasses the end of the 



linear growth phase during which electron acceleration arises (figures 2(b) and|2(c) ). The 



dominant Fourier structures of the electric field are similar to the linear phase but with 
more pronounced backward propagating harmonics. Traces corresponding to snapshots 
(Hi) and (iv) in figure 2(b) show that the electron p\\ PDF becomes increasingly flattened 



around the resonant momentum during this stage. 



From figure 3(b) we see that the peaks in energy in u, k space correspond to the 
region around ui = Qcjlh and k = ±2u pe /c. The precise values of u, k of the dominant 
waves correspond approximately to where the phase velocity of the mode matches the 
velocities of the peaks in the proton v x PDF. Thus we can infer the parallel phase 
velocity of forward v + , and backward travelling waves by reading off coordinates of 
the peaks in power from the u, k plots. Thus the electrons that are in resonance with 
these waves have forward and backward parallel momenta of 



P\l 



-v + m e / cos 9 ~ 2.5p e (t = 0) 
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p\\ - = —V-in e / cos6 ~ —3.1p e (t = 0) (11) 



in the nonrelativistic limit respectively. These momenta are indicated in figure 2(b) by 
vertical lines. These show that the instability driven by energetic protons, for these 
simulation parameters, generates waves with phase velocity resonant with the tail of the 
initial Maxwellian distribution of electron velocities, where both the gradient and the 
small number of particles are conducive to efficient Landau damping. This gives rise to 
electron acceleration and a strongly asymmetric velocity tail: hence current drive. 



5. Proton velocity space dynamics 

The particle-in-cell simulation presented here provides a full kinetic description of the 
species evolution in configuration space and velocity space. This differs from, and 
complements, analytical approaches. In the linear phase, analytical treatments are 
typically based on a small amplitude single harmonic perturbation. The linear phase 
here involves, as we have seen, fields mediating subtle interactions between the energetic 
protons, background deuterons, and electrons, in ways that would be difficult to capture 
in detail analytically. In the present section we therefore focus on the kinetics of the 
energetic proton population that drives the system, thereby obtaining a fresh perspective 
on the character of the instability during its lifetime. 

Figure [4] provides snapshots of the proton distributions in velocity space at two 
instants: (ii), during the linear phase (upper panels); and (iv), at the end of the linear 
phase (lower panels). The two left panels in figure [1] focus on velocity space. They show 
the proton velocity distribution plotted with respect to three axes defined by v» and two 
orthogonal velocity components v± t i and v± >2 in the directions perpendicular to B . Data 
plotted are for protons from a distinct region in configuration space as described below. 
The two right hand panels focus on configuration space behaviour. Each panel compares 
two subsets of protons. Each subset is localised in a distinct region of configuration space 
of extent Sx, where 5x <C A and A = 2ti / (2u pe / c) is the wavelength of the dominant 
electromagnetic wave. The two distinct regions are separated in configuration space by 
~A/2. To achieve this, protons are binned into subsets in configuration space such that 
the ratio of the size of each bin to the wavelength of the dominant electric field wave is 
Sx/\ ~ 0.057. The bins are separated by a gap of 105a; = 0.57A. 

The upper panels of figure [4] show early linear phase behaviour of the proton ring 
beam. There is a distinct wave-like structure oscillating in v» around the gyro-orbit, 
shown in panel (b) for two subsets of protons. The two datasets show velocity space 
waves that are in anti-phase with each other, consistent with resonance with a wave at 
locations A/2 apart. 
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Figure 4. Left panels: Velocities in magnetic field-aligned coordinates of protons 
found in the configuration space region < x/Ax < 16, where Ax is the spatial extent 
of one grid cell. Colour indicates speed in units ofv p (t — 0). Right panels: Proton as 
a function of gyrophase a — arctanuj^i/ux 2 f° r protons from two configuration space 
regions: < x/Ax < 16 (blue dots) and 160 < x/Ax < 176 (red pluses). Top panels 
show data from snapshot (ii), bottom panels from snapshot (iv). Velocity axes are in 



units of the proton characteristic velocity, v p (t = 0) = (l/Vp, p )£^2° 1 p Vj.Vj(£ = 0) 
where Np jP is the number of computational macroparticles representing the protons. 



1/2 



The two bottom panels of figure [4] show data from the end of the linear phase at 
snapshot (iv). The proton v\\ structure as a function of gyrophase from the same two 
configuration space bins as in panel (b) is shown for contrast in the panel (d). The 
waves in vu are beginning to break nonlinearly. 

The predominantly drift, as distinct from gyroresonant, character of the instability 
can be seen in the fully resolved velocity space of the protons. However it is challenging 
to extract the nature of the interaction between protons and the excited waves, because 
snapshots in time of proton phase space contain the history of all previous interactions. 
We now present a method for modelling the effect of previous interactions using only 
unperturbed proton trajectories which, of course, are easily calculable. 

In this linear phase the interaction between the protons and the waves is weak: 
the protons are not trapped but follow trajectories through phase space close to those 
followed by gyro-orbits in the unperturbed fields. We can therefore directly verify that 
the large scale phase space structure in figure Bib) arises from interaction with the 
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dominant wave at u ~ Qcjlh, k = —2u pe /c. To do this, we will reconstruct a phase space 
snapshot to compare with figure [4| assuming that: protons move on unperturbed gyro- 
orbits; and the perturbations in v\\ are proportional to the wave amplitude experienced 
by the proton when last at resonance. Thus we can see the history of previous wave- 
particle interactions in the phase space of protons not via perturbations to the trajectory, 
but in the history of the amplitude of the wavefield S{x, t) at the last point of resonance, 
which is calculated as follows. 

The components of parallel and perpendicular velocities of a single proton gyro- 
orbit projected onto v x can be written, referring back to ([2]), in the form 

v x (t) = ucos9 — v r (t) sin6> (12) 

Here 9 is the angle between the background magnetic field and the x-direction, and 
v r (t) = f r sin(a(t)), where the gyrophase is a(t) = Q p t + a , and a is an initial 
condition. We can calculate the time of resonance tn for this proton by noting that at 
resonance u/k = v x , so that 

v x (tR) = (jj/k = u cos 6 — v r sin(a(t R )) sin 9 (13) 

The position x(t) of the proton is given by 

x{t) = tucos9 + v r cos(a(t)) sin 9/fl p + xo (14) 

where Xq is the initial position of the proton. We eliminate the dependence on the initial 
condition xq by calculating the distance between the position of resonance xr and the 
position at time t, 

x (t) — x r — (t — tR)u cos 9 + v r sin 9(cos(a(t)) — cos(a(t#)))/fi p (15) 
The normalized amplitude of the wave at the point of resonance is 

£{x, t) = sin(fcx^(a , x, t) — ut R (a )) (16) 



We now use (15) and (16) to generate a phase space diagram for protons distributed 
across a continuum of initial gyrophase. This is shown in Figure |5j which plots the 
normalized wave field S{x^ tji) as a function of proton gyrophase a(t) at points (t 1; Xi) 
and {t\,x\ + A/2). Figure |5](a) shows how protons initially uniformly distributed in 
gyrophase resonate with the dominant backward wave at uj ~ Gujlh, k = —2u pe /c. This 
method recreates the large scale structure in Figure |4^b). The fine scale structure can 
be attributed in part to the interaction of the protons with the relatively weak forward 
travelling wave, as shown in Figure ^b) which combines the effects of both waves. Here 
&{x,t) = sin(kbXRb — UbtRb) + 0.1sin(kfXRf — 0Jft R f) is plotted at points (ti,xi) and 
{t\,X\ + A/2) against the gyrophase of the protons, where (:£&,£&) and (xf,tf) are the 
points of resonance for the backward propagating wave kb) and the weaker forward 
propagating wave (cjf, kf) respectively. In this case we take ujf/kf to equal the maximum 
in proton v x . 



Self- consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusit 




-2 2 

Gyrophase [rad] 



Figure 5. Panel (a): Normalized wave amplitude seen by protons at resonance with 
the dominant backward wave, plotted as a function of gyrophase for two regions in 
configuration space separated by half the wavelength of the dominant wave (see text). 
Panel (b): As panel (a), for the sum of the backward propagating wavemode and the 
relatively weak forward wavemode. 

6. Conclusions 

We have studied a collective instability arising from spatially localised energetic ion 
populations that can arise in tokamaks, and the phenomenology of the resultant energy 
transfer to the thermal background plasma, notably the electrons. We employ a self- 
consistent relativistic electromagnetic 1D3V particle-in-cell code to evolve minority 
energetic protons and background thermal electrons and deuterons. The 1D3V PIC 
code provides 3D velocity and electromagnetic field vectors along one spatial direction 
and time. The single spatial direction restricts the propagation direction of waves to 
an oblique angle to the magnetic field. Plasma parameters are chosen to approximate 
where possible those at the outer mid-plane edge inside the last closed flux surface in a 
large tokamak plasma. The velocity distribution of confined fusion products is modelled 
by an antiparallel travelling ring beam, motivated by observations and interpretation of 
ion cyclotron emission in JET and TFTR. 

In using a PIC code, the initial perturbation is seeded by thermal noise that 
has broad spectral content, whereas analytical studies rely on eikonal approaches. 
Computation of the proton fluctuation energy is valuable in demonstrating the existence 
of structured early time perturbations to the ensemble of particles, which are otherwise 
masked by noise. The initial perturbation to the proton distribution function generates 
a collective instability, resulting in growth of waves in which the electric field is 
substantially larger than the magnetic: predominantly electrostatic waves in the lower 
hybrid range of frequencies. Snapshots of the proton v\\ distribution function enable us 
to link the multiharmonic mode structure of the electric field to oscillations in proton 
v\\ as a function of gyrophase. 

Fast Fourier transforms of the electric field along the simulation domain before the 
linear phase, and during its early and late stages, capture proton-excited u, k structure. 



Self-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusit 



These waves drive electron acceleration through Landau damping. Directional 
asymmetry of the dominant excited wavemodes, combined with the Landau damping 
process, determines the direction of overall electron acceleration. The backward 
propagating dominant wavemode accelerates the electrons, generating a suprathermal 
tail, whereas the forward propagating wavemode interacts only weakly with the 
electrons. 

We have thus provided evidence, by direct numerical simulation, for electron current 
drive from lower hybrid waves excited by confined fusion products in a tokamak- 
relevant parameter regime. This is the first time a fully kinetic electromagnetic PIC 
code has been used to model an instability produced by velocity space population 
inversion of confined fusion products that results in electron current. These results 
provide a key building block for designing future alpha channelling scenarios, which 
may open new regimes of lower hybrid current drive where external fields and local 
fusion product population inversions together channel energy into electrons. It also 
provides an interesting instance of self organisation 51 in plasmas. 
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